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Abstract 



A new method for the computation of conserved densities of nonhnear differential- 
difference equations is apphed to Toda lattices and discretizations of the Korteweg-de 
Vries and nonlinear Schrodinger equations. The algorithm, which can be implemented 
in computer algebra languages such as Mathematica, can be used as an indicator of 
integrability. 
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1 Introduction 

Nonlinear differential- difference equations (DDEs) describe many interesting phenomena 
such as vibrations of particles in lattices, charge fluctuations in networks, Langmuir 
waves in plasmas, interactions between competing populations. Mathematically, DDEs 
also occur as spatially discrete analogues of partial differential equations (PDEs). As 
such, lattices play a key role in numerical solvers for PDEs 

In 1^, ^, ^, ^, we introduced an algorithm to find the analytical form of polynomial 
conserved densities for systems of nonlinear evolution equations. We used the concept 
of scaling symmetries or dimensional analysis. That inherently limits the algorithm to 
polynomial densities and fluxes of polynomial systems. The algorithm was implemented 
[0 in Mathematica. Here we present its extension to semi-discrete polynomial systems. 
We aim at deriving a set of independent conservation laws of DDEs, hence predicting 
integrability. 

There are several motives to find conserved densities of DDEs explicitly. The first 
few conservation laws may have a physical meaning, such as conserved momentum and 
energy. Additional ones may facilitate the study of both quantitative and qualitative 
properties of solutions 0. Furthermore, the existence of a sequence of conserved densi- 
ties predicts integrability of DDEs. Yet, the nonexistence of conserved quantities does 
not preclude integrability. Indeed, integrable DDEs could be disguised with a coordinate 
transformation so that they no longer admit conserved densities of polynomial type. 

Another compelling argument relates to the numerical solution of PDEs. In numerical 
schemes the discrete conserved quantities should remain constant. In particular, the 
conservation of a positive definite quadratic quantity may prevent the occurrence of 
nonlinear instabilities in the numerical scheme. The use of conservation laws in PDE 
solvers has been discussed in [0, |, 

For nonlinear DDEs several solution methods and integrability tests are applicable. 
The solution methods include symmetry reduction |T0|, and an extension of the spectral 



transform method |jTT|. Adaptations of the singularity confinement approach 0], the 



Wahlquist-Estabrook method []T3[, and the master symmetry technique allow one to 
test integrability of DDEs. In contrast, our method is completely algorithmic and can 
be implemented in computer algebra languages. 
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In Section 2, the algorithm is illustrated with the Toda lattice |T^. In Section 
3, we find densities for discrete analogues of the Korteweg-de Vries (KdV), nonlinear 
Schrodinger (NLS), and generalized Toda equations. These examples show subtle points 
of the algorithm. We draw conclusions in Section 4. 

2 Conserved Densities 
2.1 Definitions 

Consider a system of DDEs that are continuous in time, and discretized in the (single) 
space variable, 

Un = F(..., U„_i, U„, Un+i, ...), (1) 

where u„ and F are vector dynamical variables with any number of components. For 
simplicity of notation, the components of u„ are denoted by Un, Vn, etc. We assume that 
F is polynomial with constant coefficients. If DDEs are of second or higher order in t, 
we assume that they can be recast in the form 
For (|1|), we define a local conservation law by 

Pn = Jn ~ Jn+li (2) 

where p„ is the conserved density and J„ is the associated flux. Both functionals are 
assumed to be polynomials in u„ and its shifts. Also, (0) is satisfied on solutions of 
(|l]). Our algorithm is currently restricted to the shift-up operator [/, where (/ — f/) J„ = 
Jn — Jn+i- Minor modifications would be needed if other shift operators were used. 

Obviously, ■^{Y.n Pn) = UnPn = Hn{Jn — Jn+i)^ and this telescopic series vanishes 
if Jn is bounded for all n and J„ vanishes at the boundaries. Then, X^n Pn is constant, 
and we have a quantity that is conserved in time. 

Let D denote the shift-down operator and U the shift-up operator. Both are defined 
on the set of all monomials. If m is a monomial then Dm = m|„_»„_i and Um = 
m\n^n+i- For example, Dun+2Vn = Un+iVn~i and Uun-2Vn-i = Un-iVn- It is easy to 
verify that compositions of D and U define an equivalence relation on monomials. Simply 
stated, all shifted monomials are equivalent, e.g. Wn-itin+i = w„+2fn+4 = Un-sVn-i- 

In the algorithm below we will use the following equivalence criterion: if two mono- 
mials, nil and m2, are equivalent, mi = m2, then nii = m2 + [M„ — M„+i] for some 
polynomial M„ that depends on u„ and its shifts. For example, 'U„_2'w„ = Un-iUn+i since 

Also for later use, we call the main representative of an equivalence class, the mono- 
mial of that class with label n on m (or v). For example, UnUn+2 is the main representative 
of the class with elements m„_iu„+i, m„+im„+3, etc. We use lexicographical ordering to 
resolve confiicts. For example, UnVn+2 (not Un-2Vn) is the main representative in the 
class with elements M„_3f„_i, M„+2t'n+4, etc. 
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2.2 Algorithm 



To illustrate our algorithm, we consider the one-dimensional lattice [T^, 



jjn = exp {yn-1 - Vn) " exp (?/„ - Vn+l), (3) 

due to Toda. In (|^), y„ is the displacement from equilibrium of the rith particle with 
unit mass under an exponential decaying interaction force between nearest neighbors. 
With the change of variables, 

Un = yn, Vn = exp {yn- Vn+l), 

lattice can be written in algebraic form 

Un = Vn-l-Vn, Vn = Vn{Un - Un+l) ■ (4) 

We can compute a couple of conservation laws for (^) by hand. Indeed, iin = Pn = 
Vn-i — Vn = Jn — Jn+1 with J„ = Vn-i- We dcuotc this first pair by 

After some work, we obtain a second pair: 

n(2) = 1^ 2 , t(2) ^ 

Key to our method is the observation that (^, and @ together with the above 
densities and fluxes, are invariant under the scaling symmetry 

{t,Un,Vn) {Xt,X~^Un,X~'^Vn), (5) 

where A is an arbitrary parameter. The result of this dimensional analysis can be stated 
as follows: m„ corresponds to one derivative with respect to t; for short, Un ^ 
Similarly, v„ ~ Scaling invariance, which is a special Lie-point symmetry, is an 
intrinsic property of many integrable nonlinear PDEs and DDEs. Our algorithm exploits 
this property to find conserved densities, which now proceeds in three steps. 

Step 1: Determine the weights of variables 

The weight, w, of a variable is by definition equal to the number of derivatives with 
respect to t the variable carries. Weights are positive, rational, and independent of n. 
We set w{^) = 1. In view of (|^), we have w{un) = 1, and w{vn) = 2. 

The rank of a monomial is defined as the total weight of the monomial, again in 
terms of derivatives with respect to t. For instance, the rank of each monomial in p^^^ 
is two. Observe that in each equation of (^, all the terms (monomials) have the same 
rank. This property is called uniformity in rank. Densities and fluxes are also uniform 
in rank, and from (|^), it follows that rank(J„) = rank(p„) + 1, since vj{-^) = 1. 
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Conversely, requiring uniformity in rank for each equation in (Q) allows one to com- 
pute the weights of the dependent variables. Indeed, w{un) + 1 = w{vn),w{vn) + 1 = 
w{un) + w{vn), yields w{un) = 1, w{vn) = 2, which is consistent with (|]). 

Step 2: Construct the form of the density 

As an example, let us compute the form of the density of rank 3. List all monomials in 
Un and Vn of rank 3 or less: Q = {Un^, UnVn, Un, Vn}- 

Next, for each monomial in Q, introduce enough t-derivatives, so that each term 
exactly has weight 3. Thus, using (H), 



d° 

— (u 3- 
dt° ^ " ' 


) = Mn^ 


d° 




d^("") 








d2 




d/ N 

1 = ^(^n-l - Vn) = Un- 





Gather the resulting terms in a set = {ur? ,UnVn-i,UnVn,Un-iVn-i,Un+iVn\ ■ Identify 
members that belong to the same equivalence classes and replace them by the main 
representatives. For example, since UnVn-i = Un+iVn both are replaced by UnVn-i- 
Doing so, Ti is replaced by X = {u^ , UnVn-i, UnVn}, which contains the building blocks 
of the density. Linear combination of the monomials in X with constant coefficients q 
gives the form of the density: 

Pn = CiUn^ + C2UnVn-l + CsUnVn- (6) 

Step 3: Determine the unknown coefficients in the density 

Now we determine the coefficients Ci through C3 by requiring that (0) holds. During this 
step we also compute the unknown flux J„. 

Compute pn using (^. Then use (^) to remove ?)„, etc. After grouping the terms 

Pn = (3Ci - C2)Un^Vn-l + (Cs - 3Ci)Un^Vn + (Cs - C2)f„-lU„ 

Use the equivalence criterion to modify p„. For instance, replace Un-iUnVn-i by UnUn+iVn+ 
[un-iUnVn~i — WnMn+i'^^n] • The goal is to iutroduce the main representatives. Therefore, 

Pn = (3Ci - C2)Un^Vn-l + (C3 - 3Ci)'U„^i)„ 

+ (C3 - C2)VnVn+l + [(C3 - C2)Vn-lVn - (C3 - C2)VnVn+l\ 
+C2UnUn+lVn + [c2M„_lM„t;„_l - C2U„M„+it;„] 
+C2Vn'^ + \c2Vn-\ " C2Vn \ - C^UnUn+lVn " CsVn^ . 
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Next, group the terms outside of the square brackets and move the pairs inside the 
square brackets to the bottom. Rearrange the latter terms so that they match the 
pattern [J„ — J„+i]. Hence, 

Pn = (3Ci - C2)Un^Vn-l + (C3 - 3Ci)Un^Vn 

+ (C3 - C2)VnVn+l + (C2 - C3)m„M„+iV„ + (C2 - C3,)Vn^ 
+ [{(C3 - C2)Vn^iVn + C2U„_iM„t;„_i + C2i'„_l^} 
-{(Cs - C2)VnVn+l + C2UnUn+lVn + C2vJ}]. 

The terms inside the square brackets determine: 

Jn = (cs - C2)Vn-lVn + C2M„_iti„t;„_i + C2W„_l^. (7) 

The terms outside the square brackets must all vanish, yielding 

S = {3ci - C2 = 0, C3 - 3ci = 0, C2 - C3 = 0}. (8) 

The solution is 3ci = C2 = C3. Since densities can only be determined up to a multi- 
plicative constant, we choose ci = |, C2 = C3 = 1, and substitute this into and (0). 
Hence, 

Pn = I Un + Un{Vn-l + Vn) , Jn = Un-lUnVn-l + Vn-l^ ■ 

Analogously, we computed conserved densities of rank < 5 for (^). They are: 

Pn^ = Un, p'i^ = \Un + Vn, P^f^ = \Un + Un{Vn-l + Vn), 

Pn^ = lUn^ + Un'^{Vn-l + Vn) +UnUn+lVn + J + VnVn+l, 

p^^^ = + uJ{Vn-l + Vn) + Mn^^n+l^n(^n + '"n+l) 

+UnVn~l{Vn~2 + Vn~l + Vn) + M„W„(w„_i + Vn + Vn+l) ■ 

Ignoring irrelevant shifts in n, these densities agree with the results in [1T6| . 
To illustrate how our algorithm works for DDEs with parameters, consider 

Un = aVn-l-Vn, Vn = Vn {P - Un+l) , (9) 

where a and f3 are nonzero parameters. In it was shown that (|^) is completely 
integrable if a = /5 = 1. 

Using our algorithm, one can easily compute the compatibility conditions for a and 
(3, so that (P) admits a polynomial conserved densities of, say, rank 3. The steps are the 
same as for (^). However, (|^) must be replaced by 

S = {3aci — C2 = 0, /3c3 — 3ci = 0, ac^ — C2 = 0, /?C2 — C3 = 0, ac2 — C3 = 0}. 

A non-trivial solution 3ci = 02 = 03 will exist if and only if a = f3 = 1. 

Analogously, (|) has density p^^^ = of rank 1 if a = 1, and density p^^^ = + 
of rank 2 if a /3 = 1. Only when a = f3 = 1 will have conserved densities of rank > 3. 
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3 Examples 

3.1 The Volterra Equation 

Consider the integrable discretization of the KdV equation: 

tin = Un {Un+l - Un-l) , (10) 

which is known as the Kac-Van Moerbeke equation or a special form of the Volterra 
system. It arises in the study of Langmuir oscillations in plasmas, and in population 
dynamics fT^, |TB|, |g. 

Notice that ([l0|) is invariant under the scaling symmetry [t, Un) —>■ (At, A~^m„). Hence, 
Un corresponds to one derivative with respect to t, i.e. Un ^ All terms in ([10|) have 
the same rank if w{un) + 1 = 2w{un), thus, w{un) = 1, which agrees with the scaling 
symmetry. 

Let us find the form of density with rank 3. Forming all monomials of m„ with rank 3 
or less yields the list Q = {un^ , tt„^, Un}. Introducing the necessary t-derivatives, leads 
to H: 

r 3 2 2 2 2 1 

\Un ,Un Un+li^n-lUn ,tf„M„_(_i , W^-ilt^li^+i ,M„M„_|_itt„_|_2 ,tt„_2li„_ili„, tt„|. 

Using Un-2Un-lUn = Un-lUnUn+1 = M„'Un+l'Un+2 , ^in-l^n^ = li„U„+l^ and ti„_i^-U„ = 

Un^Un+1, we obtain the list 

7" _ r 3 2 2 -I 

— \Un ,Un Un+l,UnUn+l , UnUn+lUn+2 f ■ 

A linear combination of the terms in X with constant coefficients Cj gives 

Pn = Cl Un^ + C2 UnUn+1 + C3 M„M„+1^ + C4 UnUn+lUn+2- 

Proceed with step 3. After differentiation, shifting and regrouping 

Pn = (3Ci - C2)ujun+1 + (C3 - 3Ci)^l„M„+i^ + 2(c2 - Cs)Un^Un+l'^ 
+2(C3 - C2)UnUn+l^Un+2 + (C2 " C4)^Z„^U„+iU„+2 

+ (C4 - C3)m„M„+iM„+2^ + [Jn " Jn+l], (H) 

with 

Jn = -{3CiUn-lUn^ + 2C2^X„-lU„^M„+i + CsUn-lUnUn+l^ + C4M„_lM„U„+iM„+2) • 

The monomials outside the square brackets in (0) must vanish. This yields 

S = {3ci - C2 = 0, C3 - 3ci = 0, C2 - C3 = 0, C2 - C4 = 0, C4 - C3 = 0}. 
Choosing ci = |, one has C2 = C3 = C4 = 1. Therefore, 

Pn = I Un + UnUn+1 [Un + Un+1 + Un+2) , 
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Analogously, for ( [T0| ) we computed the densities of rank < 5 : 



Pn^ = i^^*^ + + l^n^^n+l^ + Mnttn+l^(^n+l + Wn+2) 

[Un + Wn+i + M„+2) + UnUn+lUn+2^{Un+2 + Un+s) 
+UnUn+lUn+2Un+3{Un + ^n+l + + + Mn+4)- 

3.2 Discretizations of the nonlinear Schrodinger equation 



In P0| , , Ablowitz and Ladik studied properties of the following integrable discretiza- 
tion of the NLS equation: 

± U*^Un{Un+l + Un-l) , (12) 

where is the complex conjugate of m„. We continue with the + sign; the other case 
is analogous. Instead of splitting m„ into its real and imaginary parts, we treat m„ and 
= as independent variables and augment ([T2|) with its complex conjugate equation. 
Absorbing i in the scale on t, we get 

Un = Un+l - '2Un + Un-1 + UnVn{Un+l + ^n-l), 

Vn = -{Vn+1 - 2Vn + Vn-l) - M„W„(w„+i + (13) 

Since Vn = u*^, we have wivn) = w{un)- 

Neither of the equations in (|I3p is uniform in rank. To circumvent this problem we 
introduce an auxiliary parameter a with weight, and replace ( pTSD by 

Un = a{Un+l-2Un + Un-l) +UnVn{Un+l + Un-l), 

Vn = -a{Vn+l - 2Vn + Vn-l) " UnVn{Vn+l + Vn-l). (14) 

Uniformity in rank requires that 



w{Un) + 1 = w{a) + w{Un) = 2w{Un) + Uj{yn) = 3w{u 
w{Vn) + 1 = W{a) + w{Vn) = 2w{Vn) + w{Un) = 3w{Vn), 



which yields w{un) = w{vn) = \,w{a) = 1, or, Un^ ~ ~ a ~ ^. 

Recall that the uniformity in rank requirement is essential for the first two steps of 
the algorithm. However, after Step 2, we may set a = 1. The computations now proceed 
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as in the previous examples. We list some conserved densities of (|T3|), which correspond 
with those in [pO[] : 



pI^^ = CiUnVn-l + C2UnVn+l, 

Pn^ = CiilUn^Vn-l^ + UnUn+lVn-lVn + UnVn-2) 

+ C2{\UnVn+l^ + M„M„+1 V„+l V„+2 + UnVn+2), 

pf^ = Ci[lUn^Vn-l^ + UnUn+lVn-lVn{UnVn-l + Un+lVn + Un+2Vn+l) 

+ UnVn-l{UnVn-2 + Un+lVn-l) + M„t^„ (m„+1^^„-2 + Un+2Vn~l) + ^n^^n-s] 

+ C2[lujvn+1^ + UnUn+lVn+lVn+2{UnVn+l + 'U„+iW„+2 + U„+2i'„+3) 

+ UnVn+2{UnVn+l + Un+lVn+2) + (^n+l^^n+l + Un+2Vn+2) + UnVn+s]- 



As shown in scheme {^), if defined on an infinite interval, admits infinitely many 
independent conserved densities. Although it is a constant of motion, we cannot derive 
the Hamiltonian of (p!2D, for it has a logarithmic term 

We also computed conserved densities of the non-integrable standard second-order 
scheme [0, 

i iln = Un+1 - 2Un + + S'U*^^, (15) 

for the NLS equation. Instead of (|^) one has 

Vn = -aivn+i - 2vn + v^-i) - 2m„z;^. 

Here, ~ ~ a ~ ^. We could only find two independent conserved densities. 
Indeed, after setting a = 1, 

pll^ = UnVn, P^n^ = UnVn + UnVn-l + Mn^n+l- 

3.3 Generalized Toda lattices 

Recently, Suris [^, showed the integrability of the chain 

Vn = 2/n+ie^^"+^"^"^ - e'^^y"+'-y"^ - y„_ie(^"~^"-i) + e^^y"-y"-'\ (16) 

which is related to the relativistic Toda lattice. With the change of variables, Un = 
i/n, Vn = exp {un+i - Vn) , (HD cau be written as 

iln = Vn{Un+l - Vn) - Vn-l{Un-l - Vn^l) , Vn = Vn{Un+l - Un) ■ (17) 

Here, Un ^ Vn ^ ^, and we computed five conserved densities for (|1^). The first three 
are: 

Pn ^ ^™ ^™ 1 Pn ^ i 

Pn^ = ki^n^ + 2^^n^) " '"n(^^n-l^ + ^n^) + UnUn+lVn- 
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Suris 1^ also investigated 

ijn = Vn [exp {yn+i - Vn) " exp {yn - Vn-i)] , (18) 

which is closely related to the classical Toda lattice (|^). The same change of variables 
as for (plGf) allows one to write (ffl) as 



iln = Un{Vn - Vn~l) , Vn = Vn{Un+l - Un) ■ (19) 

Again, Un ~ ~ and the first four conserved densities are 

P^^^ = Un + Vn, p^^'^ = ^{Un^ + Vn^) + Un{Vn-l + Vn) , 

Pn^ = liUn^ + Vn^) +Un'^iVn-l+Vn) +UniVn-l'^ +Vn'^) +UnVniVn-l + Un+l), 
Pn^ = liuj + Vn'^) + uJ{Vn-l + Vn) + Un{Vn-l^ + vj) + lUn^{Vn-l^ + Vn'^) 

+UnUn+lVn{Un + ^n+l) + 2M„Vn («n^^n-l + Un+lVn) 

+«„i)„_iZ;„(i)„_i + Vn) + UnUn+lVn{Vn-l + Vn+l) ■ 

4 Conclusions 

We developed a Mathematica program, called diffdens.m, and used it to compute all 
the conserved densities presented in this paper. For lattices with parameters, the code 
automatically determines the compatibility conditions on these parameters so that a 
sequence of conserved densities might exist. 

The existence of a large number of conservation laws is an indicator of integrability 
of the system. Therefore, by generating the compatibility conditions, one can analyze 
classes of parameterized DDEs and filter out the candidates for complete integrability. 

Future generalizations of the algorithm will exploit other symmetries in the hope to 
find conserved densities of non-polynomial form. 
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